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Abstract 

Keratin are among the most abundant proteins in epithelial cells. Functions of the 
keratin network in cells are shaped by their dynamical organization. Using a collection 
of experimentally-driven mathematical models, different hypotheses for the turnover 
and transport of the keratin material in epithelial cells are tested. The interplay 
between turnover and transport and their effects on the keratin organization in cells are 
hence investigated by combining mathematical modeling and experimental data. 
Amongst the collection of mathematical models considered, a best model strongly 
supported by experimental data is identified. Fundamental to this approach is the fact 
that optimal parameter values associated with the best fit for each model are 
established. The best candidate among the best fits is characterized by the disassembly 
of the assembled keratin material in the perinuclear region and an active transport of 
the assembled keratin. Our study shows that an active transport of the assembled 
keratin is required to explain the experimentally observed keratin organization. 


Introduction 


The epithelial cytoskeleton is characterized by abundant keratin intermediate filaments 
(Fig.[T]). The cytoplasmic keratin filament network is responsible for the mechanical 
stress resistance of epithelial cells and contributes significantly to epithelial 


stiffness 1^ 28 . The importance of keratins for epithelial tissue stability is reflected by 


a large group of genetic skin blistering diseases that are caused by point mutations in 
keratin-encoding genes [^[^. The mechanical functions not only rely on static 
resilience but also necessitate a high degree of plasticity, for example in migrating cells 
during wound healing [^. The current view is that keratins act as general stress 
absorbers protecting epithelial cells not only against mechanical insults but also against 
irradiation or osmotic and microbial challenges. Thus, keratins are involved in heat 
shock response, apoptosis and organelle homeostasis [^. Furthermore, functions 
affecting processes such as proliferation, differentiation and inflammation are also 


dependent on keratins (see recent reviews in 11 ^). 
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Figure 1. Keratin network. Image taken from a time-lapse fluorescence recording of 
a single hepatocellular carcinoma-derived PLC cell of clone PK18-5 


30 stably 


expressing fluorescent fusion protein HK18-YFP consisting of human keratin 18 and 
enhanced yellow fluorescent protein. Bar 10 \xm,. 


All of these functions are tightly coupled to the keratin network dynamics (see |S1 


Video). Examination of cultured epithelial cells producing fluorescent keratins has 


provided evidence for the different mechanisms that are involved in the continuous 
renewal and reshaping of the keratin system 


16 . On the basis of these observations, a 


biological model of the keratin cycle has been proposed by Leube and Windoffer et al. 


17 33 . This biological model takes into account the assembly/disassembly and 


transport of keratins. In this study, it was proposed that assembly and disassembly 
occur in topologically defined regions with assembly taking place predominantly in the 
cell periphery while disassembly takes place primarily in the perinuclear region. The 
biological model further postulates active transport of insoluble assembly stages of 
keratins toward the nucleus and rapid diffusion of soluble subunits throughout the 
cytoplasm. To the best of our knowledge, how these different processes are coupled and 
regulated is not yet known. 

In a previous work we developed methods to examine and quantify the keratin 
transport and turnover in epithelial cells 21 ; the spatial distribution of the assembled 


keratin material in epithelial cells was available at 24 hours and 48 hours after seeding. 
The current effort is to use the available qualitative and quantitative observations to 
derive, from first principles, experimentally-driven mathematical models that could 
yield hypothetical predictions testable in laboratories. Our approach is unique, it 
translates experimental observations and data into a series of alternative plausible 
mathematical models or scenarios to further advance our understanding of the critical 
parameters in keratin cycling. Fundamental to this approach is the fact that optimal 
parameter values for each scenario are established and out of this set, a single scenario 
is identified that best fits experimental observations and data. 

Hence, a collection of mathematical models resulting from different assumptions of 
the keratin transport, assembly and disassembly is designed to investigate the effects of 
the interplay between turnover and transport on the keratin organization. The 
collection is built as a well-designed scientific experiment by considering control and 
knockout of processes. To highlight and confirm the importance or existence of a given 
process, scenarios in which the process is absent are also considered and tested. Model 


29 














responses are then compared to experimental observations and data published in 21 


to 


identify optimal parameter values that yield the best fit of each of the models to the 
experimental data. Finally, we identify, using an information-theoretic approach, the 
best scenario or model given the data and candidate models under study. 

By employing this reductionist phenomenological approach and systematic 
evaluation of the different scenarios we not only confirm the proposed transport features 
of the keratin cycle and the restricted disassembly in the perinuclear region but also 
find that the assembly throughout the cytoplasm fit best to the experimental data. 
Furthermore, our particular approach allows us to demonstrate that the inward motion 
experimentally observed is not an emergent behavior but is an inherent property of the 
keratin material organization and it is due to an active transport thereby confirming 
recent experimental observations 


17 33 


Methods 


Experimental data 

In Moch et al. , the spatial distribution of the assembled keratin material in 
epithelial cells is measured for 15 minutes at 24 hours and 48 hours after seeding. The 
shape of each epithelial cell is normalized to fit a circle of fixed radius . The average 
normalized spatial distribution is calculated over 50 cells at 24 hours (resp. 84 cells at 
48 hours). The average speed and direction of the motion of the assembled keratin 
material are measured and determined at every location within the normalized cell. 
Finally, at every spatial location, the net assembly/disassembly is calculated. Hence, 
regions with preferential assembly and disassembly are identified. We will refer to 
regions of assembly as Sources and regions of disassembly as Sinks. More details on the 


experimental data can be found in 10 21 


In the present work, cells are represented as one-dimensional cross-section domains. 
A diameter of the normalized cell is used as the spatial domain which is centered at the 
center of the cell and is of length 2L with L = 22.5/im. Moch et al. recorded the 
fluorescence intensity of fluorescent protein-labeled keratins in cells. Assuming a 
proportionality between fluorescences and concentrations, and knowing from the 
mean concentration for keratin in keratinocytes, we convert fluorescence intensities to 
concentrations {^M) as follows: Concentration = Fluorescence x (Mean Concentration 
/ Mean Fluorescence). In Fig. the average spatial distribution, the speed of the 
assembled keratin, regions of assembly (Sources) and regions of disassembly (Sinks) on 
the one-dimensional cross-section domain are displayed. 


Mathematical models 

To study its organization in cells, the keratin material is categorized into a soluble 
pool composed of the soluble keratin, and an insoluble pool representing the 
assembled keratin material. The state variables used to represent the soluble and 
insoluble pools are: 

• S{x,t) denotes the concentration of the soluble keratin material at position x at 
time t, 

• I{x,t) denotes the concentration of the assembled keratin material at position x 
at time t. 

From here onwards, the one-dimensional spatial domain representing the cell is defined 

by 

H = {x : —L < X < L}, 
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Figure 2. Experimental data adapted from |2l| . 2.1: The spatial distribution of 
the assembled keratin material. The mean concentration of the keratin material used in 
the conversion of fluorescence intensities to concentrations is estimated to be equal to 
520/iM in [^. Circles represent raw experimental data; curves are fits of experimental 
data. Details on P{x) and ffinai{x) can be found in Appendix 4 2.2: The speed of the 
assembled keratin material at 24 hours. 2.3: Regions of assembly and disassembly 
denoted Sources and Sinks respectively. In all figures, the cell is represented by a 
one-dimensional cross-section domain centered at the center of the cell; plasma 
membrane positions are at ±L = ±22.5^m, the nuclear envelope is located at ±7.5^m 
and the center of the cell is located at zero. 


where x = 0 is the center of the cell and x = ±L are the boundary positions of the 
plasma membrane (Table [^. The general framework of the model, derived from first 
principles based on experimental observations, takes into account the turnover and 
transport for both soluble and insoluble pools and can be stated verbally as: 

{ Rate of change of 5" = Transport of S' — Assembly of S in / -|- Disassembly of I in S, 
Rate of change of / = Transport of / -I- Assembly of S in / — Disassembly of I in S. 

Hence, the generalized model has the following expression; stated mathematically as: 

{ f)Q 

— = Ts{S)-A{S)+V{I), 

( 1 ) 

TiiI)+A{S)-ViI), 

where Ts{-) (resp. T/(-)) is the functional term describing the transport of the soluble 

d / 

pool (resp. insoluble pool). For example, TsiS) = — — Js{x,t) (resp. 

T/(/) = — — J/(x,t)^ where Js{x,t) ^resp. Ji{x,t)'^ describes the flux of the soluble 

pool (resp. insoluble pool) at position x from the left (x = —L) to the right (x = L) at 
time t. The function A(S') (resp. D(/)) is the assembly term (resp. disassembly term). 
To investigate the interplay between turnover and transport on the organization of the 
keratin material several assumptions are proposed for each of the functionals T 5 , Tj, A 
and V. 
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Table 1. Model parameters 


Parameter Definition 


Value (Unit) 


Reference 


L 


fo{x) 


Half-length of spatial do- 22.5 (/j.to) 
main H 

Initial distribution of as- Eq. 0 if^M) 
sembled keratin in cells at 
24 hours 


21 


21 


Ds 


Dr 


}[x) 


u(x) 


Diffusion coefficient of sol- 0.88 ± 0.08 {[imr/s) 

uble pool 

Diffusion coefficient of in- < 10~^Ds /s) 

soluble pool 

Speed of insoluble pool 0.002 — 0.008 {^m/s) 


16 


Almost constant speed of Eq.(12) (fim/s) 
active transport of insolu¬ 
ble pool 


21 


34 


Variable speed of active Eq.(13) (fim/s) 
transport of insoluble pool 


Rate of assembly of soluble TBD (linear model: s 

pool nonlinear model: 


21 


^ass (^) 


ks 


kdis 


kdis (:^) 


ki 


Localized rate of assembly Eq.(14) with kass TBD 

(linear model: nonlin- 


of soluble pool 


ear model: ^M.s 


Concentration for half- TBD (/rM) 
saturation of assembly rate 
(nonlinear model) 

Rate of disassembly of in- TBD (linear model: 
soluble pool into soluble nonlinear model: 
pool 


Localized rate of disassem- Eq.(15) or ( |16| ) with kdis 
bly of insoluble into soluble TBD (linear model: 
pool nonlinear model: 


Concentration for half- TBD (pM) 
saturation for disassembly 
rate (nonlinear model) 

The model parameters used in all the scenarios. (TBD = To Be Determined by fitting 
model solutions to experimental data). 


Modes of transport Molecules of the soluble pool are assumed to be subjected to 
the Brownian motion; the soluble pool is diffusible with a diffusion coefficient Ds- The 
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functional term of transport for the soluble pool is given by 


T,(S) := D.0, 


( 2 ) 


Passive transport for both pools is assumed in all models. Diffusion is assumed for the 
insoluble pool to describe the wiggling motion of the keratin filaments in cells. The 
diffusion coefficient of the insoluble pool is set to be much smaller than that of the 
soluble pool: 0 < Dj < 10~^Ds. It is assumed that only the insoluble pool can be 
driven by an active transport. Experimental evidence show that the assembled 
intermediate filament proteins in the form of filament precursors, squiggles or filaments 
move along microtubules and actin filaments by interacting via molecular 
motors [^[^|^[^|^. An active transport (also called the inward drift) of the 
insoluble pool from the plasma membrane to the center of the cell is hypothesized based 
on reports of the motion of the assembled keratin material mostly towards the nucleus 
in epithelial cells 15 33 34 . The speed of the active transport is set to be almost 
constant v{x) everywhere or variable u{x) everywhere (Fig. [^. Based on experimental 
observations, both speeds are assumed to decay around the nucleus towards the center 
of the cell. The variable speed u{x) is estimated using the profile of the average speeds 
measured in (Fig. [^2). The magnitude of the almost constant speed v(x) is set to 
be the average value of the variable speed over the cell. Details on the derivation of the 
estimates of v{x) and u{x) are given in Appendix 1 Hence, the functional term for the 


transport of the insoluble pool can take three different forms: 


n ^ 

Tiil):={ Di^ 


no drift, 


92 / dl 

7 :-^ + sgn{x)v{x)^—, inward drift with almost constant speed, (3) 
ox‘^ ox 

d'^I dl 

D /—-77 + sgn{x)u{x)^:—, inward drift with variable speed, 

aa;2 

where v(x) is given in (12) and u(x) in (13); the two functions representing the speeds 
of the active transport are graphed in Fig. ^ The function sgn{x) defined by 


sgn{x) = 


1 a: > 0 , 




describes the inward direction of the active transport at any location of the spatial 
domain D centered at zero. 

Combining the modes of transport for the soluble and insoluble pools, three types of 
transport are considered for the keratin material in the epithelial cell. 


Expressions of assembly / disassembly reactions In the present work, the 
turnover is composed of two reactions: the assembly / aggregation / polymerization of 
units of the soluble pool to grow the insoluble pool and the disassembly / solubilization 
/ depolymerization of the insoluble pool into units of soluble pool. It is assumed that 
the assembly process is a function of the soluble pool only, whereas the disassembly 
process is a function of the insoluble pool only. The simplest case is to consider a linear 
model that assumes linear exchanges between the two pools. The linear model can be 
stated as follows: 

±(A{S)-V{I)) ■.= ±(kass{-)S-kdrs{-)l), (4) 
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Location (p, m) 

Figure 3. Profiles of the active transport speeds u{x) and v{x) for the 
insoluble pool both of which are derived from the experimental measurements in 21 


The function u{x) is derived from the profile of speeds and v{x) is approximated as the 
average value of these speeds. Details on the derivation of u(x) and v{x) are given in 
Appendix 1 


or 


where kass{‘) (resp. kdis{-)) is the rate of assembly of the soluble pool (resp. 
disassembly of the insoluble pool). Both rates can either be constant, kass and kdi. 
space-dependent kass{x) and kdis{x). 

In the second case, the turnover is assumed to depend on the enzymatic 
activities [^. For instance, the solubilization of the insoluble pool into soluble proteins 
(disassembly) is triggered by a kinase activity and the assembly of insoluble pool is 
induced by the dephosphorylation of soluble proteins by a phosphatase 
turnover term is assumed to be of the Michaelis-Menten form stated as: 


12 . The 




= ± 


^ass {■)S kd^s{■)I 
ks + S ki + I 


( 5 ) 


where kassi’) (resp. kdisi’)) is the maximum rate of assembly of the soluble pool (resp. 
disassembly of insoluble pool) and ks (resp. kj) is the concentration at which the 
assembly (resp. disassembly) rate is half of kass{-) (resp. kdis{-))- When 
Michealis-Menten dynamics are used, the model is called nonlinear. Both rates can 
either be constant or space-dependent describing the intracellular localization of the 
post-translational modification enzymes. 


Profiles of assembly and disassembly rates As previously mentioned, the 
assembly and disassembly rates, fcoss(’) and kdisi’), used in the linear and nonlinear 
models can be constant or space-dependent functions. The profile (shape) of the 
space-dependent function kassix) is derived from the spatial profile of regions of 
assembly (Sources) measured in m (Fig- §3). Details of the derivation of the Sources 
kassix) are given in Appendix 2| 

Two types of shapes for kdisix) are assumed to represent two types of localization of 
the disassembly in cells. First, similarly to the assembly rate, the profile of the 
disassembly rate is deduced from the experimental data published in (Fig. [^3); the 
spatial profile of the disassembly regions. Sinks, is used to build the shape of the first 
space-dependent disassembly rate. This disassembly rate is called of type Sinks. Second, 
it is assumed that disassembly is localized around the nucleus; a mollified step-function 
is designed to describe this assumption. This second space-dependent disassembly rate 
is called of type Mollify. Details on the derivation of the two kdis ix) rates are given in 
[Appendix 3[ 
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For the sake of illustration, the two profiles of kass{') and the three profiles of kdisi') 
are shown in Fig. 
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Figure 4. Possible profiles for the assembly and disassembly rates. 4.1: 

kass = 10”^ when the linear model is used (resp. when the nonlinear model 

is used) and kass{x) is computed using (14|. 4.2: kdis = 10”^ when the linear model 
is used (resp. when the nonlinear model is used) and kdis{x) are computed 

using ( |T^ for Sinks and (16) for Mollify, respectively. Details on the derivation of 
knasix) and kdisix) rates are given in Appendix 2 and [Appendix 3[ respectively. 


Parameter values in (14), (15) and (16) are chosen in such a way that their profiles give 


the same total amount of assembly / disassembly over the spatial domain D. 


Accounting for the three modes of transport and considering the turnover described 
by either linear and nonlinear models with 6 possible combinations of the profiles for the 
assembly and disassembly rates, a collection of 36 scenarios (models) is defined; each 
scenario follows the general form stated by system ([^. Details of the 36 scenarios are 
given in Fig. All scenarios are considered with the same initial conditions given by 

fs'(a:,to) =7^Mx), 

< for all x G ft, (6) 

I/(a;, to) = foix), 


where fo{x), defined in (0, is chosen to be a mollified version of the profile of the 
assembled keratin at time to = 24 hours averaged over all the normalized cells (Fig. |^. 
Details of the derivation of fo{x) are given in [Appendix 4 Initial conditions describe 
the observed fact that the soluble pool (resp. insoluble pool) represents 5% (resp. 95%) 
of the total keratin material [^. All scenarios are considered with boundary conditions 
describing the impermeability of the plasma membrane for the keratin material 


= Q = t > to, (7) 

where Js is the flux of the soluble pool and Jj is the flux of the insoluble pool as defined 
below system (0. The model parameters used in all the scenarios are listed in Table 0 
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None 



Diffusion — with 

constant 
speed 


with 

variabie 

speed 



Linear 


Nonlinear 


Linear 


Nonlinear 


Linear 


Nonlinear 





Passive Active Type of 
transport transport exchange 


Non-localized ^ 

Localized - Source ^ 
Non-iocalized ' 
Localized - Source • 
Non-localized ' 
Localized - Source 

Non-localized ^ 

Localized - Source > 
Non-localized 

Localized - Source ■ 
Non-localized ' 

Localized - Source 

Assembly 





Non-localized (1) 
Localized - Sinki 
Localized 


- Siiikjf2) 

- Mollify { 


3) 


!^o°c"aiPzl"c]'^®^i^ft5) 

Localized - Mollify (6) 
Non-localized (7) 

tsiiiiiSra/®’ 

Non-localizec 
Localized - S 
Localized - Mollify (12) 

Non-localized (13) 
Localized - Sink n41 
Localized - Mollify (15) 

Localized - Mollify (18) 

Non-localized (19) 
Localized - SinkizO) 
Localized - Mollify (21) 

Localized - Mollify (24) 
Non-localized (25) 

- Localized - Sihk f26) 

• Localized - Mollify (27) 

Non-locajizei 


Localized - SIhk (29) 
Localized - Mollify (30) 

Non-localized (31) 
Localized - Sinkj[32) 
Localized - Mollify (33) 

Localized - Mollify (36) 

Disassembly 


Scenario i 

Diffusion 

No drift Cst drift Var. drift 

Lin. Nonlin. 

Cst Ass. Sources 

Cst Diss. Sinks Mollify 

1 

1 

1 

1 

1 

1 

2 

1 

1 

1 

1 

1 

3 

1 

1 

1 

1 

1 

4 

1 

1 

1 

1 

1 

5 

1 

1 

1 

1 

1 

6 

1 

1 

1 

1 

1 

7 

1 

1 

1 

1 

1 

8 

1 

1 

1 

1 

1 

9 

1 

1 

1 

1 

1 

10 

1 

1 

1 

1 

1 

11 

1 

1 

1 

1 

1 

12 

1 

1 

1 

1 

1 

13 

1 

1 

1 

1 

1 

14 

1 

1 

1 

1 

1 

15 

1 

1 

1 

1 

1 

16 

1 

1 

1 

1 

1 

17 

1 

1 

1 

1 

1 

18 

1 

1 

1 

1 

1 

19 

1 

1 

1 

1 

1 

20 

1 

1 

1 

1 

1 

21 

1 

1 

1 

1 

1 

22 

1 

1 

1 

1 

1 

23 

1 

1 

1 

1 

1 

24 

1 

1 

1 

1 

1 

25 

1 

1 

1 

1 

1 

26 

1 

1 

1 

1 

1 

27 

1 

1 

1 

1 

1 

28 

1 

1 

1 

1 

1 

29 

1 

1 

1 

1 

1 

30 

1 

1 

1 

1 

1 

31 

1 

1 

1 

1 

1 

32 

1 

1 

1 

1 

1 

33 

1 

1 

1 

1 

1 

34 

1 

1 

1 

1 

1 

35 

1 

1 

1 

1 

1 

36 

1 

1 

1 

1 

1 


Figure 5. The 36 scenarios to be considered. Top: Number in parentheses is the 
scenario index i. Bottom: The numerical value 1 denotes that the process of interest is 
in Scenario i. 


Comparison between mathematical models and experimental 
data 


Parameter estimation Let p denote the set of all the model parameters for each 
scenario. To estimate the optimal set of parameter values p for each scenario, the 
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Figure 6. Initial profile of the insoluble pool fo{x) defined in 

solution of each scenario is compared to experimental data using the following objective 
function (an estimation of the distance between the experimental data and the model 
response): 


Wx 

^(p) — ^ ^ 5 t/ina/ 5 p) 


( 8 ) 


where tfinal is equal to 48 hours. The experimental data at tfinal is represented by 
ffinal {x) and the model solution at tfinal of the considered scenario evaluated with the 
parameter values p is I{x,tfinahP)- To obtain the model solution for each scenario, the 
corresponding system is numerically integrated from to = 24 hours to 50 hours using the 


MATLAB solver for partial differential equations, pdepe 20 . Solutions are computed at 
Nx locations of the spatial domain with = 200. To carry out computations, the 
raw data in which the concentration of the assembled keratin material is known at 623 


spatial locations is approximated by ffinai{x) defined in Appendix 4 by (18); ffinai{x) 
is the fit of the average profile of the assembled keratin measured after 48 hours on the 
normalized cells (see the black profile in Fig. [^1). 

The estimation of the parameter values for the 36 scenarios, i.e. the determination of 
parameter values p that provide the best fit to the experimental data, is done by 
minimizing the objective function, 4'(p), such that 


4' (p) = min 4/ (p). 

p 


( 9 ) 


Minimization of the objective function is done by a parallelizable genetic algorithm 
described in [^. 

In this study, only constant parameter values of the turnover reactions are optimized. 
When the linear model is considered, for each scenario, two parameters kass and kdis 
are estimated. When the nonlinear model is considered, four parameters, kass, ks, kdis 
and fc/, are estimated for each scenario. In all scenarios, the diffusion coefficients for the 
soluble and insoluble pools as well as the active transport speeds v{x) and u{x) are 
considered as fixed/known parameters or functions and are set to values measured in 
previous studies [^21 ^ (Table [^. The diffusion coefficients are taken as 
Ds = 0.88pm^/s and Dj = 9.5 x 10~'^Ds, and the “constant speed” u in v{x) as 


defined in (12) is set to be u = 0.0025pm/s 34 . 


Model selection In order to select the best model out of the 36 scenarios, we use an 
information-theoretic approach. Specifically, we employ the Akaike information criterion 
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(AIC) [ 3,13 to select from all the scenarios, the best model that best captures 
experimental observations given the collection of models considered in this study. Since 
we use the Least Squares principle to estimate the values of the constant free 
parameters (2 free parameters for linear models and 4 for nonlinear models), the Akaike 
criterion for Scenario i, AlC'i, takes the following form: 


AlCi = In (a^) + 2K, (10) 

where is the estimate of the variance with \E'i(p) the residual estimated 

for Scenario i and the size of the sample (A(e=200). K is the number of estimated 
parameters; i.e. the number of free parameters in Scenario i plus one for the estimate of 
the variance {K = 3 for linear model and AT = 5 for nonlinear models). The scenario 
with the lowest AIC value is the best model. The AIC selects a model with the least 
number of parameters that best fits experimental observations. To rank and compare 
scenarios the Akaike weights Wi are calculated and these are known as the weights of 
evidence in favor of Scenario i being the actual best model given the experimental data 
and the collection of scenarios considered. The Akaike weights are expressed as follows: 


exp(-l/2A,) 

Ef=iea;p(-l/2Ar)’ 


with Ai = AlCi — min AlCi 

i 


( 11 ) 


where R represents the number of scenarios considered (i? = 36) and A^ being the 
difference in AIC with respect to the AIC of the preferred scenario min AlCi = AlCmin- 

i 

It must be noted that the Akaike weights Wi sum to I and are interpreted as the 
probability that Scenario i is the best model given the experimental data and the 
collection of scenarios considered. The models, ranked from the largest to the smallest 
Akaike weights, whose Akaike weights sum to 0.95 form the confidence set of the models 
that captures, more faithfully, experimental data. The ratio of the Akaike weights 
Wi/wj (also known as the evidence ratio) is used to compare model pairs. Furthermore, 
the relative importance of a process can be estimated by summing the Akaike weights of 
all scenarios involving the process of interest. We will denote by w'^ the sum of the 
Akaike weights of the scenarios including the process of type *. This sum can be 
interpreted as the probability that the process of type * is the best type of the process 
given the experimental data and the collection of models considered. 


Results 


Best 


scenario 


We employ a two-step process in order to find the best scenario. First, the best fit for 
each scenario is found by minimizing the objective function Secondly, the Akaike 
information criterion (101 is used to select the best of the best fits out of all the 
scenarios. 

The best fit for each of the 36 scenarios is presented in Fig. in the following order: 

1. Row-wise: 


• Scenarios in the first three rows (1 to 3) have linear terms for turnover. 

• Scenarios in the last three rows (4 to 6) have Michaelis-Menten type turnover 
terms. 

• Rows 1 and 4 are scenarios with constant disassembly rate. 

• Rows 2 and 5 are scenarios with localized disassembly rate of type Sinks. 


11/29 









• Rows 3 and 6 have scenarios with localized disassembly rate of type Mollify. 

2. Column-wise: 

• The first two columns (1 and 2) include scenarios with diffusion alone 
without drift. 

• Columns 3 and 4 display scenarios with drift with an almost constant speed. 

• The last 2 columns (5 and 6 ) are scenarios with drift with variable speed. 

• Odd columns (1, 3 and 5) are scenarios with constant assembly rate. 

• Even columns (2, 4 and 6 ) are scenarios with localized assembly rate of type 
Sources. 

According to the AlCi values, the best model, the best of the best fits, is identified 
as being Scenario 21 , AIC 21 = min (3’’'^ column of Tablej^. Since Scenario 21’s 

Akaike weight is over 0.95 (5*^ column of Table [^, it is the only model out of the set 
considered that satisfies the confidence criterion; Scenario 21 matches very well 
experimental observations and data. The second best model is Scenario 31. The 
evidence ratio of scenarios 21 and 31 W 21 /UI 31 is equal to 36.67 x 10®; i.e. Scenario 21 is 
about 36 millions times more likely than Scenario 31 to be the best model given the 
experimental data and the collection of models considered. We consider this to be 
strong evidence in support of Scenario 21. 
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Table 2. Results of the model selection for the best fit of each of the 36 
scenarios. 


Scenario i 

K 

AIC, 

A, 

Wi 

Rank 

1 

3 

2012.996 

510.3519 

0 

13 

2 

3 

2019.144 

516.4999 

0 

19 

3 

3 

2018.634 

515.9895 

0 

18 

4 

3 

2512.249 

1009.604 

0 

36 

5 

3 

2283.252 

780.6081 

0 

28 

6 

3 

2340.069 

837.4254 

0 

31 

7 

5 

2011.436 

508.7923 

0 

11 

8 

5 

2011.724 

509.0804 

0 

12 

9 

5 

2013.609 

510.9647 

0 

14 

10 

5 

2014.842 

512.1975 

0 

15 

11 

5 

2015.148 

512.5043 

0 

16 

12 

5 

2017.034 

514.3897 

0 

17 

13 

3 

1839.930 

337.2861 

0 

7 

14 

3 

2479.026 

976.3817 

0 

35 

15 

3 

1738.738 

236.0939 

0 

5 

16 

3 

2114.919 

612.2750 

0 

23 

17 

3 

2458.497 

955.8532 

0 

34 

18 

3 

2002.972 

500.3279 

0 

9 

19 

5 

1668.771 

166.1274 

0 

3 

20 

5 

2387.721 

885.0770 

0 

32 

21 

5 

1502.644 

0 

0.99999 

1 

22 

5 

2073.809 

571.1653 

0 

21 

23 

5 

2338.356 

835.7119 

0 

29 

24 

5 

1893.875 

391.2309 

0 

8 

25 

3 

1703.773 

201.1292 

0 

4 

26 

3 

2410.565 

907.9213 

0 

33 

27 

3 

2006.784 

504.1403 

0 

10 

28 

3 

2067.488 

564.8442 

0 

20 

29 

3 

2243.066 

740.4220 

0 

27 

30 

3 

2145.639 

642.9953 

0 

26 

31 

5 

1537.479 

34.83477 

2.7 X 10"“ 

2 

32 

5 

2339.200 

836.5561 

0 

30 

33 

5 

2119.447 

616.8033 

0 

24 

34 

5 

2074.025 

571.3811 

0 

22 

35 

5 

1834.186 

331.5420 

0 

6 

36 

5 

2143.595 

640.9505 

0 

25 


The numerical value 0 in the 5*^ column denotes a value lower than 10“^^. In the 6*^ 
column, “Rank” denotes the ranking of scenarios, i.e. the descending order of Akaike 
weights Wi- The use of weights allows the ’’quantitative” comparison of the adequacy of 
scenarios because of the normalization to 1. 


Scenario 21 is characterized by an inward drift with an almost constant speed for the 
insoluble pool, turnover terms of Michaelis-Menten type, a constant assembly rate and a 
disassembly rate of type Mollify. The profiles of assembly and disassembly rates are 
displayed in Fig. For this scenario, the estimated optimal parameter values are 
fcass = 9.3819/rM/s, ks = 570.73^M, k^ig = 0.9998/iM/s leading to kmax = l-976^M/s 
in (16) and fc/ = 976.07. Based on the Michaelis-Menten constants for the assembly 
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and disassembly processes, since ks < kj, an enzyme that would be involved in the 
solubilization of the assembled keratin material requires a higher substrate 
concentration to achieve a given reaction speed than an enzyme that would be involved 
in the assembly of soluble proteins. Snapshots of the soluble and insoluble pool profiles 
taken every 30 minutes from Iq to tfinal are displayed in Fig. After only 2 hours, the 
solution stabilizes to its final profile. Scenario 21 preserves the repartition of the keratin 
material between the soluble and insoluble pools over time, about 95% of the keratin 
material is assembled to form the insoluble pool. Finally, the characteristic time scales 
of the passive transport {T^iffusion), active transport (rDrift) and turnover (tR eaction) 
for Scenario 21 are estimated by using an adimensionalization o f system ([T|) 
corresponding to Scenario 21. Details on calculations are given in [Appendix 5| The time 
scales of the processes included in Scenario 21 are ordered as follows: 

TReactioni^^ s) < '^Diffusioni^^ '®) ^ ^Drifti^O s) < s). 

Using Peclet’s number {Pe = Tuiffusion/Torift ^ 1), it is found that the transport of 
the assembled keratin material by drift is faster than by diffusion, the dominant mode 
of transport is the active transport. Using Damkohler number 

{Da = TRrift/TReaction ^ !)> it is found that the active transport time scale is greater 
than the reactive time scale; overall. Scenario 21 is controlled by active transport. 

The complete ranking based on the Akaike weights uji of all the scenarios is given in 
the 6*^ column of Table It is worthwhile remarking that differences in A/C^, A^, 
could have been used for ranking purposes; in this case, we would have obtained the 
same ranking. The rationale of using uji is that it allows us to quantify how preferably 
each candidate model is via the normalization to 1. 

Importance of the type of process 

Using the sums of the Akaike weights, , we investigate several questions to evaluate 
the relative importance of the different types considered for a process (Tables and |^. 
The set of the 36 scenarios is partitioned into categories with respect to the types of 
processes considered. For instance, considering the transport process, the collection of 
the 36 scenarios is partitioned into 3 categories (Fig. [^: scenarios with no drift 
(Scenarios from 1 to 12), scenarios with an almost constant speed drift (Scenarios from 
13 to 24) and scenarios with a variable speed drift (Scenarios from 25 to 36). The sum 
of the Akaike weights of each category is then calculated, compared and ordered to 
characterize which type is more likely to be present. In what follows the sign > denotes 
the relative importance, measured in probabilistic terms, of the type of process. 


Table 3. Importance of the type of process. 


Type of Transport 

Type of Turnover Term 

Type of Assembly 

Type of Disassembly 

^ NoDrift 
=0 

^CstDrift 
= 1 

^Var Drift 

= 2.7 X 10"^ 

^Linear 
= 0 

^ N onlinear 

= 1 

^CstAss 
= 1 

^Source 

=0 

'^CstDis 

= 2.72 X 10"® 

Ws^nk 

=0 

"^tloll 
= 1 


denotes the sum of the Akaike weights of scenarios including the type * of the process. Type of transport 
- No drift for the insoluble pool (No drift) vs Drift with almost constant speed (Cst drift) vs Drift with 


variable speed (Var. drift). Type of turnover terms - Linear vs Nonlinear. Type of assembly - Non-localized 
(constant) vs Localized of type Sources. Type of disassembly - Non-localized (constant) vs Localized of type 
Sinks vs Localized of type Mollify. 
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Table 4. Importance of the type of process. 


Scenario i 

Cst Ass. Dis. 

Cst Ass. Sinks 

Cst Ass. Moll. 

Sources Cst Dis. 

Sources Sinks 

Sources Moll. 

1 

1 






2 


1 





3 



1 




4 




1 



5 





1 


6 






1 

7 

1 






8 


1 





9 



1 




10 




1 



11 





1 


12 






1 

13 

1 






14 


1 





15 



1 




16 




1 



17 





1 


18 






1 

19 

1 






20 


1 





21 



1 




22 




1 



23 





1 


24 






1 

25 

1 






26 


1 





27 



1 




28 




1 



29 





1 


30 






1 

31 

1 






32 


1 





33 



1 




34 




1 



35 





1 


36 






1 


stAssDis 

= 2.7 X 10"® 

C stAssSink 

=0 

'^CstAssMoll 
= 1 

^ SourceDis 
=0 

wi c ■ 1 

bourcebink 

=0 

Source. Ad oil 

=0 


Combinations of assembly/disassembly rate profiles - Constant assembly and disassembly rates vs constant 


assembly rate and disassembly rate of type Sinks vs constant assembly rate and disassembly rate of type 
Mollify vs assembly rate of type Sources and constant disassembly rate vs assembly rate of type Sources and 
disassembly rate of type Sinks vs assembly rate of type Sources and disassembly rate of type Mollify. Top: 
The numerical value 1 denotes that the combination of processes of interest is in Scenario i. Bottom: w'^ 
denotes the sum of Akaike weights of scenarios including the process combination as indicated at the top. 


• What is the most likely type of transport for the keratin material in 
cells given the experimental data and the model collection considered? 

Is it that there is no drift of the insoluble pool (diffusion only for both pools), drift 
with an almost constant speed for the insoluble pool or drift with a variable speed 
for the insoluble pooU Each category includes 12 scenarios (see 3’’'^ to 5*^ columns 
in Fig. [^. From the Akaike weights, it is obvious that scenarios including 
diffusion only (with no drift) are not supported by experimental data since 
w'^oDrift ~ ^ (Table 1^. Diffusion alone is not enough to explain the experimental 
data. Given the data and the model collection, the type of transport can be 
ordered as follows (Table [^: 

Drift with almost constant speed > Drift with variable speed ^ No drift. 


15/29 














• What is the most likely type of turnover between the soluble and 
insoluble pools given the experimental data and the model collection 
considered? Is linear exchange or Michaelis-Menten type (underlying an enzyme 
activity) more likely! Each category includes 18 scenarios (see 6*^ to 7*^ columns 
in Fig. [^. From Tableenzymatic activities are more likely to be present. Hence 

Nonlinear exchange ^ Linear exchange. 

• What is the most likely rate profile of the assembly process of the 

keratin material in cells given the experimental data and the model 
collection considered? Is a constant assembly rate or assembly mainly localized 
at the cell membrane periphery more likely! Each category includes 18 scenarios 
(see 8*^ to 9*^ columns in Fig. [^. From Table the non-localized rate of assembly 
is preferred to the rate of the localized assembly; . Hence 

Non-localized assembly 3> Assembly localized at cell membrane periphery. 

• What is the most likely rate profile of the disassembly process of the 
keratin material in cells given the experimental data and the model 
collection considered? Is a constant disassembly, disassembly of type Sinks or 
disassembly localized around the nucleus more likely! Each cate gor y is composed 
of 12 scenarios (see 10*^ to 12*^ columns in Fig. [^. From Table]^ none of the 
scenarios that include the disassembly rate of type Sinks is supported by 
experimental data. 

Disassembly localized around the nucleus > Non-localized disassembly ^ 

Sinks-type profile. 

• What is the most likely combination of the assembly and disassembly 
rate profiles given the experimental data and the model collection 
considered? Is a combination of constant assembly and disassembly rates, a 
constant assembly rate and disassembly rate of type Sinks, a constant assembly 
rate and disassembly rate of type Mollify, an assembly rate of type Sources and 
constant disassembly rate, an assembly rate of type Sources and disassembly rate 
of type Sinks or an assembly rate of type Sources and disassembly rate of type 
Mollify more likely? Now, interactions between the assembly and disassembly 
processes are investigated. Six categories are defined. In each category, there are 6 
scenarios as listed in the top of Table Overall, as in Scenario 21, the 
combination of non-localized assembly and disassembly localized around the 
nucleus is preferred (Table |^. 

Non-localized assembly and disassembly localized around the nucleus > 
Non-localized assembly and disassembly ^ Any other combination. 


Discussion 


The primary aim of the present work is to investigate, through mathematical modeling 
driven by experimental observations, mechanisms contributing to the organization of 
the keratin material in epithelial cells and subsequently to test the biological model for 
the keratin dynamics proposed by Leube and Windoffer et al. in 2011 


17 33 . From 


first principles, we formulate a collection of mathematical models capturing various 
combinations of biological processes describing the spatio-temporal dynamics of the 
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keratin material in epithelial cells. By using techniques in parameter estimation, we find 
optimal reaction kinetic parameter values for each model that best captures 
experimental observations. We go one step further and employ an information-theoretic 
approach for model selection to determine the best model among all the best fit models 
that captures the key processes of the experimental data. 

As previously highlighted, the model framework and the description of processes 
considered are driven by experimental data and observations. Experimental data used in 
this study provide the spatial distribution of the assembled keratin concentration at 24 
and 48 hours. As only macroscopic information on the keratin organization is available, 
the models only describe the keratin material as soluble or assembled and consider 
exchanges between these soluble and insoluble pools combined with transport events. 
For instance, as experimental data identify the existence of regions with preferential 
assembly and disassembly (Fig. [^3), the assembly and/or disassembly processes are 
assumed to be localized in some scenarios. When observing the perpetual inward motion 
of the keratin network (see SI Video| , the “natural assumption” for the transport of the 
assembled keratin would be the existence of an inward active transport. However, we 
decide to consider scenarios with only passive transport (no active transport). Why? 
First, it is well known that combining appropriate reaction terms and diffusion (only) 
can lead to the emergence of complex behavior such as traveling waves and/or pattern 
formation (see some examples in |^). Secondly, to reinforce the predictive power of our 
study. In our approach, we do not only design mathematical models but also calibrate 
models (parameter estimation) and then evaluate (model selection) how each model 
performs. To highlight and confirm the importance or existence of a given process, 
scenarios in which the process is absent must also be considered and evaluated. Our 
strategy of considering “unrealistic scenarios” and “realistic scenarios” in the collection 
of models to be evaluated mimics the biological experimental protocols such as knockout 
and controls. For instance, it has been shown that the keratin assembly / disassembly 
depend on post-translational modifications of keratins due to enzymatic activities 12 


In our study, a total of 36 scenarios that divide into 18 scenarios with a turnover of type 
linear (no enzymatic activities) and 18 scenarios with a turnover of type nonlinear 
(enzymatic activities) are investigated. There is a one-to-one correspondence between 
the 18 scenarios of the 2 groups. When model performances (Akaike weights) are 
compared, scenarios having a nonlinear turnover perform better than the corresponding 
ones with a linear turnover. Hence, the enzymatic activity is detected by the model 
selection as existing but also preferable. This outcome allows us to judge that our 
approach combining the mathematical modeling and model selection performs correctly, 
it gives us the confidence and trust in our conclusions. 

Following this methodology, it turns out that Scenario 21 was the best of the best 
fits. Scenario 21 is in good agreement with the biological model proposed in 


17 33 


Scenario 21 and the biological model share the following common key features: diffusion 
of the soluble pool, inward active transport of the insoluble pool and disassembly of the 
insoluble pool localized around the nucleus. Moreover, as experimentally observed. 
Scenario 21 preserves the repartition of the keratin material over time; the keratin 
material is mainly insoluble in epithelial cells. Both the proposed biological model and 
Scenario 21 hypothesize an inward active transport for the assembled keratin material. 
ExaminiM the Akaike weights of the models with no active transport (Scenarios 1 to 12 
in Table0 I'** and 2"“^ columns in Fig. [^and I'** to 3’^'^ columns in Table|^, we go one 
step further and show that the active transport is a requirement to explain the 
experimental data and that the experimentally observed inward motion of the 
assembled keratin is not an emergent phenomenon but it is due to an active transport. 
Furthermore, comparing the characteristic time scales of processes, we establish that the 
keratin dynamics are mainly controlled by the active transport in Scenario 21. 
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Interestingly, Scenario 21 concurs with the proposed biological model about the 
existence of a disassembly localized in the perinuclear region. An experimental protocol 
must now be developed to further work out the details of this conclusion. It is worth 
pointing out that some characteristics of the proposed biological model were not tested 
due to the form of the mathematical models. For instance, the scenarios considered here 
describe the turnover in terms of the assembly of soluble proteins and disassembly of 
aggregated proteins. In the biological model, the nucleation of filaments, i.e. the 
initiation of filaments, is explicitly described and localized at the cell periphery. 
However, the nucleation process is not explicitly described in any of the present models. 
It is worth mentioning that at the beginning of this work a nucleation term was 
included in the models; the effect of this term did not change significantly the results; 
hence, to reduce the complexity of models and the number of parameters the nucleation 
term was subsequently dropped. Investigation, by mathematical modeling, of the 
nucleation process is currently being carried out in a separate and on-going study. 

When only scenarios with a variable drift are considered (Scenarios 25 to 36, 5*^ and 
6 *^ columns in Fig. [^, the preferred combination of assembly and disassembly rate 
profiles is a constant rate of assembly and disassembly as in Scenario 31, which is the 
second best model. The best model, Scenario 21, and the second best model, Scenario 
31, share the non-compartmentalization of the assembly process. However, in Scenario 
21 the active transport has an almost constant speed whereas in Scenario 31 the speed 
is variable. For the disassembly, this process is non-localized in Scenario 31 whereas it is 
localized at the perinuclear region in Scenario 21. The spatial variability of the speed 
u{x) in Scenario 31, in particular, the almost-zero speed at the nucleus locations (see for 
X G [—7.5, 7.5] in Fig. “compensates” for the non-compartmentalization of the 
disassembly process. After comparing the features of Scenarios 21 and 31, we inspect 
their profiles. The profile obtained with Scenario 21 fits very well the experimental data 
on non-perinuclear locations (see for x ^ [—7.5, 7.5] in Fig. [^33). On the other hand, 
the profile resulting from Scenario 31 fits very well experimental data on perinuclear 
locations (see for x G [—7.5, 7.5] in Fig. [^23). As an improvement of Scenario 21, we 
would expect that a wider decay in speeds around the nucleus modeled, for instance, 
with a smaller value of a in (12) would result in a better fit of the perinuclear region. 

Scenarios 29 and 35 include the variable speed measured from experimental data, a 
localized assembly rate having the spatial profile deduced from Sources and a localized 
disassembly rate whose shape is derived from the profile of the Sinks. All the 
characteristics extracted from the experimental data are included in Scenario 29 (linear 
model) and Scenario 35 (nonlinear model). If one wanted the best model that takes into 
account all the experimental features, then Scenario 29 or 35 would be the best scenario. 
Scenario 29 is ranked 27*^ and Scenario 35 is ranked 6*^ (Table [^. Similarly to the 
general trend, for the same set of assumptions, the use of the Michaelis-Menten type 
turnover term generally provides a better representation of experimental data than the 
use of the linear turnover term. The evidence ratio of Scenarios 35 and 29 uj 35 /uj 2 g is 
ridiculously large; Scenario 35 is much more adequate to represent the experimental 
data than Scenario 29. However, Scenario 35 is still only ranked 6*^. The 
failure/mismatch of Scenario 35 might be explained by the redundancy of the 
information existing in the variable speed and the net assembly/disassembly region 
profiles extracted from the experimental data. Furthermore, according to the Akaike 
weights, disassembly rates of type Sinks (the type deduced from experimental data) are 
less likely to occur given the experimental data and the collection of models considered. 
A similar conclusion is obtained for the assembly rate of type Sources. This could be a 
consequence of our too conservative interpretation of the regions of preferential 
assembly or disassembly. In our assumptions, in zones of preferential assembly, the 
disassembly rate is set to be about null and vice versa (Fig. §3). 
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In summary, to model the keratin dynamics in epithelial cells we characterized the 
keratin material into two pools, the soluble and the insoluble pool; and, the events 
considered are the turnover and transport for both pools. The modeling assumptions 
used, for instance, the diffusion of the soluble pool, are based on biological observations 
and experimental data. The collection of the models considered in this study is designed 
to answer a set of questions such as “what is the mode of transport of the keratin 
material in epithelial cells?”. After optimizing parameter values, model selection and 
evaluation methods applicable to non-nested models are used to discriminate between 
the candidate models, identify the best model and quantify how models under 
consideration are adequate to explain the experimental data. Note that the ranking of 
the models (scenarios) and the relative importance of the different types of processes 
(Tables [2]|^ are only valid in the context of the experimental data and the set of the 
candidate models considered here. For instance, considering other hypotheses such as 
the non-negligence of the anterograde motion of the assembled keratin material or the 
stabilization (protection against disassembly) of the keratin filaments involved in the 
nuclear cage would have led to a different collection of scenarios in which our best 
scenario could have failed to be the best one. Furthermore, as some modeling 
assumptions are directly derived from experimental data, missing information in 
experimental data might have been prejudicial for the correct approximation / modeling 
of processes; for instance, missing information about the speed of assembled keratin in 
the cell periphery result in small values for the variable speed close to the plasma 
membrane (Fig. [^2). Keeping in mind the limitations of our approach, important 
conclusions have been reached such as 


an active transport of the assembled keratin material is required, thereby 


confirming recent experimental observations 17,33 


enzymatic activities regulating the assembly / disassembly are more likely to 
occur, 

the assembly process is more likely to be non-compartmentalized in cells, 

last but not least, a unique best model strongly supported by experimental data is 
identified, this scenario is in good agreement with the biological model previously 
proposed in 


17 33 . Interestingly, the best scenario supports the perinuclear 


localization of the disassembly process hypothesized in the biological model. 


Mathematical models of the keratin intermediate filament organization in cells were 
previously proposed [^[^[^[^|^; however, none of these models included the effects 
of transport of the assembled material on its organization. More importantly, only the 
behavior of models were validated qnalitatively, no comparisons to experimental data 
were carried out. It is worth noting that other studies of the intermediate filaments 
dynamics combining mathematical modeling and experimental data approaches were 
carried out but on neurofilaments in neurons (see for example [^|18|). 


Conclusion 


Given the experimental data published in 21 , through modeling and simulations, we 


investigate the effects of the interplay between turnover and transport on the keratin 
spatio-temporal organization in epithelial cells. Out of all the scenarios investigated, a 
scenario strongly supported by experimental data is found that best captures most of 
the hallmarks of the experimental observations. This scenario predicts the diffusion of 
soluble keratin, an inward active transport of the assembled keratin and the disassembly 
localized around the nucleus triggered by enzymatic activities as well as the assembly 
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process that is non-compartmentalized over the cell. The value of our models is 
reflected in their predictive nature; first, the approach predicts the localized disassembly 
at the perinuclear region and second, that the experimentally observed inward motion is 
not an emergent behavior but that it is an inherent property of the organization of the 
keratin material in epithelial cells and it is due to an active transport. 


Appendix 1 

Based on experimental observations, the speed of the assembled keratin material is 
assumed to decay to almost zero around the nucleus. The nuclear envelope positions are 
at a; = ±7.5^m. In order to describe the decay of the speed around the nucleus such 
that it is almost constant, the following function is used: 

v(x) = u(l — exp(—ax^)), (12) 


with a = 0.05 and u being in the range given in Table [l] For numerical simulations, u is 


set to 0.0025/j,m/s which represents the average value of speeds measured in 21 


For the variable speed case, a symmetrical function over the spatial domain 12 is 
sought to describe the speed u(x). Averaging over the symmetrical spatial locations 
values of the average speed measured in (Fig. [^2) and curve fitting with a sum of 
Gaussian functions of these values, an estimate of u(x) is obtained and expressed as 
follows: 


2 

u(x) = '^ai X exp(-((x - hi)/cif), (13) 

i=l 


with the following coefficients: oi = 0.003372, hi = 17.39, ci = 7.577, 02 = 0.003378, 

62 = —17.41 and ci = 7.546. The estimate u{x) is almost zero at the nucleus-locations 
and is symmetrical around the center of the cell (Fig.[^. 

For numerical simulations, when needed, the following function 


s{x) 


2 

1 -I- e-^/“ 


- 1 


with a = 1 is used as a “smooth analogue” of sgn{x). 


Appendix 2 

To obtain the space-dependent function kass (x) the profile of the regions of assembly 
(Sources) published in (Fig. [^3) is first made symmetrical by averaging values at 
the symmetrical spatial locations and then the symmetrical Sources profile is fitted 
using the sum of Gaussian functions: 

3 

kassix') — kmax ^ ^ ^ CXp(^ ((x h^^ j ) F ^6aseZme- (^'^) 

i=l 


The best fit is obtained with the following coefficients: ai = 1.684 x 10°, bi = 20.13, 

Cl = 2.105, 02 = 8.493 x 10^ 62 = -19.04, ca = 1.374, 03 = 1.333 x 10^ 63 = -20.74, 
and C 3 = 1.734. The value of the coefficient kmax determines the maximal value of the 

^I'i.kass ^baseline) 


peaks. When kmax = 


with k, 


baseline 


= X 10 the total 


V^J2i=iai\ci\ 

amount of assembly over the cell is the same as that of the case of a constant assembly 
of level kass- The value kass in kmax is determined by fitting model solutions to 
experimental data. An illustration of the shape of kassix) is given in Fig. 10 
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Appendix 3 

To obtain the first space-dependent function kdis{x), the profile of regions of 
disassembly (Sinks) published in (Fig. [^3) is made symmetrical and fitted using the 
sum of Gaussian functions: 


kdis{^^ —kn 


8 

2=1 


Oi X exp{-{{x - hi)/cif). 


(15) 


The best fit of the symmetrical profile of Sinks is obtained with the following 
coefficients: oi = -3.383 x 10^ bi = -11.51, ci = 4.158, as = -1.004 x 10®, 

62 = -14.03, C 2 = 2.236, 03 = 2.445 x 10®, 63 = -11.39, C 3 = 4.339, 04 = -2.032 x 10^, 

64 = 5.922, C 4 = 5.385, 05 = -2.058 x 10^ 65 = -6.077, C 5 = 5.778, ag = 3.793 x 10®, 

bQ = 11.37, C 6 = 4.322, 07 = -1.008 x 10®, bj = 14.03, C 7 = 2.238, ag = -4.742 x 10®, 

6g = 11.46, and Cg = 4.196. The coefficient k^ax determines the amplitude of the peaks. 

1! d/‘^ s 

When kmax = g , the total amount of disassembly over the cell is the same 

V^T,i=la^\ci\ 

as that of the case of the constant function of level kdis- The value kdis in kmax is 
determined by fitting model solutions to experimental data. An illustration of the shape 


of kdis{x) obtained from Sinks is given in Fig. 11 


A second function is hypothesized to represent the localized disassembly kdis{x) 
around the nucleus. A mollified piecewise function is used and is of the form 


kdisi^^') — 


kbas eline ? 

h 

^m.n 


x^- 


4e(ai - 02 ) 

— (e — ai)^h 


kmax{^ ai) 
- X 

2e(ai — 02 ) 

,X 4:C{Qj\ G,2)k}iaselir\ 


^baseline 

h 

^m.n 


4e(ai - 02 ) 
u a: - ai 

‘^max ? 

a2 - ai 

2 kmaxi^ 0 , 2 ) 

rX + —— - ^X 


4e(a2 - Oi) 2e(a2 - Oi) 

((e - 02 ^ + 4eai)fc max + 4e(ai - a 2 )haseii 


4e(a2 - ai) 

^max ^ 3 ) 


^baseline 4“ ^ 

h 

^max 

46(03-04)“' ' 26(03-04)“" 

^ ((e + 03 )^ — 4604)^ max -I- 46(o3 - ai)kbaseli 

46(03 - 04 ) 

kh r +k — k ^ ~ 

^baseline ' ^max ^max ? 

04 — 03 

^max 2 kyriax\^ O 4 ) 


46(04 - as) 

(e + a 4 )‘^k„ 


26(04 - as) 

X -f 46 (o4 as)kbaseline 


khaselinet 


4e(a4 - as) 


X < ai — e, 

ax — e < X < ai + €, 
ai + e < X < a 2 — £, 

02 — 6 < a; < 02 -b 6, 
02 -b 6 < a: < 03 — 6, 

03 — 6 < a; < 03 -b 6, 
03 -b 6 < a: < 04 — 6, 


04 — 6 < a; < 04 -b 6 , 
04 -b 6 > a:. 


(16) 


where oi = —15/rm, 02 = —7.5pm, as = 7.5pm, 04 = 15pm and 6 = 1. When 
kmax = 2{kdis - haseHne) with kbaseiine = kdis X 10 “^, the total amount of disassembly 
over the cell is the same as that of the case of a constant disassembly of level kdis', kdis 
is determined by fitting the model solutions to experimental data. An illustration of the 
shape of (16) is given in Fig. 
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Appendix 4 

To define the initial condition fo{x), the experimental profile of the assembled 
keratin material measured at 24 hours is used (gray circles in Fig. il)- The polynomial 
8 

P{x) = ^^pix' for which pg = 5.441 x 10“®, pj = 3.397 x 10“^^, pe = —5.379 x 10“®, 
1=0 

P 5 = -2.077 X 10"^®, p 4 = 0.01062, pg = 2.801 x lO'^® p 2 = 0.4104, pi = 1.4 x 10"^^ 
and po = 506.5 fits well the data at 24 hours (gray curve in Fig. [^1). The function P{x) 
is then modified to obtain a function satisfying the boundary conditions defined in 0. 
The initial condition /o(x) describing the profile of the insoluble pool at 24 hours is 
expressed as follows: 


fo{x) = < 


P(ai+e)-c 2 , (e - ai)(P(ai + e) - c) 

-7^5-^ ”1-XT-^ 

4e^ 2e^ 

(e — ai)^P(ai + e) + ( 3 e^ + 2 aie — af)c 


+ 


P{x), 


4^2 


a: < oi — e, 

ai — e < a; < oi + e, 
tti + e < a; < 02 — e. 


P(a2 — e) — c + 2eD{a2 — e) 2 
-TT- ^ 

4g2 

(e — a2)(P(a2 — e) — c) — 2ea2D{a2 — e) 


(17) 


+ 


+ 


2e2 

( 3 e^ + 2026 — a2)P(a2 — e) 

4g2 

2e(e^ — a2)D{a2 — e) + (e — 02)^0 

4^2 > 


c, 


02 — e < a; < 02 + e, 
02 + e < a;, 

8 


with oi = —21fj,m, 02 = —oi, e = 1, c = 50pM and D{x) = ’'^^ipiX^ The function 


i=l 


/o(x) is graphed in Fig. 

The average profile of the assembled keratin material computed on the normalized 


cells after 48 hours of seeding 21 is represented by: 

4 


ffinal ^ 


(18) 


i=0 


where p 4 = -0.003255, pg = 2.61 x lO'^^ p 2 = 0.4899, pi = 1.558 x lO'^® and 
Po = 604.1. ffinai{x) is graphed in Fig.j^l (black curve). 


Appendix 5 

The system corresponding to Scenario 21 is now nondimensionalized to estimate the 
characteristic time scale of each of the processes involved. Consider new independent 
variables z and r and new dependent variables S{z,t) and I{z,t) defined such that: 

X = Az, t = Bt, S{x,t) = sS{z,t), I(xA) = 

Scenario 21 takes then the following form: 


dS^ 

dr 

m 


BDs d^S 



^dis (4z)/ 
ki/i + I 


.S 


BcDsd'^i Bu , A \di B ( 


42 


kdrs{Az)i 
kiji + J 
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with f{Az) being the functional component of v{x) defined as 
f{Az) = (1 — exp{—a{Az)‘^)) with a = 0.05 and kdis{Az) being the function kdis{-) 
defined in (16) and evaluated at Az. Taking A = £ (£ = 2L), B = £/ m , s = ks and 
i = ks then 

£ 

X = £z and t = —r 
u 

and the adimensional version of the system is 


55 1 d^S 

dr ePe dz"^ 


— Da 


5 kdis^Az) / k ass^ 


1 + 5 


ki/ks 


91 ^ 9'^i , , 5 / ^ 


-I 

5 


kdis{Az)/k 

ass^ 


1 + 5 ki/ks + I 


is the Peclet number of the insoluble 


The constant Pe = 

TDrift £/u eDs 

pool that compares active transport and diffusion processes. The characteristic length £ 
is chosen to be equal to the length of the spatial domain fl, £ = 2L. The Damkohler 


number Da = 


'^Drift 


compares reaction and active transport processes. The 


'^Reaction 

characteristic active transport time is Thrift = £lu as defined in the Peclet number. 

k 

The characteristic reaction time is Tueaction = 1/k for which n = -. Note that 

ks 

sgn{z)f{Az) is about ±1. 

Recall the parameter values obtained from the parameter estimation: 
kass = 9.3819^M/s, ks = 570.73/iM, kj = 976.07/j,M/s and kmax = 1.976/iM/s in (16) 
that correspond to a constant disassembly rate of level kdis = 0.9998/iM/s. Hence, from 
the definition of kmaxi kdis{Az) can be approximated by kdis = 0.9998. The values of 
the fixed parameters are Ds = 0.88iJ,m^/s, Dj = eDs with e = 9.5 x 10“"^ and 
u = 0.0Q25fim/s. The following estimates are then obtained: 


• Adimensional parameters: kdisiAz)/kass = kdis/kass = 0.1 and kj/ks = 1.7; 

• Diffusion time scale: Toiffusion = e/D. 

- For the soluble pool: Xoiffusion = ^l^s ~ 2.3 x lO^s, 

- For the insoluble pool: ~ 2.4 x 10®s; 

• Active transport time scale: Thrift = £lu = 1.8 x lO^s (time the fluid needs to 
flow through the characteristic length); 


• Reaction time scale: ks/kass ~ 61s (time for reaction to equilibrate). 


To determine what mode of transport (passive vs active) and what process 
(transport vs reaction) dominate the dynamics of the keratin material, the Peclet and 
Damkohler numbers are used: 


• Pe= {u£)/{eDs) ~ 135 S> 1. Diffusion is small compared to active transport for 
the insoluble pool, the transport of the assembled keratin material by drift is 
faster than by diffusion, active transport is the dominant mode of transport; 

• Da = {£kass)/{uks) « 296 ;+ 1. Active transport time scale is greater than the 
reaction time scale; the overall process in controlled by active transport. 
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Supporting Information Legends 


SI Video 

Dynamics of the keratin network in a cell. Time-lapse fluorescence microscopy of 
hepatocellular carcinoma-derived PLC clone PK18-5 stably expressing fluorescent fusion 


protein HK18-YFP 30 depicting the dynamic properties of the keratin filaments over a 


time period of 15 hours. Bar 10 /rm. 
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7.31 Scenario 9 7.32 Scenario 12 7.33 Scenario 21 7.34 Scenario 24 7.35 Scenario 33 7.36 Scenario 36 

Figure 7. Best fit for each of the 36 scenarios. Red curve is the response Ii{x,tfinal,p) at 48 hours of Scenario i. Black circles are data at 
48 hours as in Fig. [^1. 
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Figure 8. Details about the best model Scenario 21. 8.1: Snapshots from 24 
hours (gray) to 48 hours (black) taken every 30 minutes of the profile of the soluble pool 
over the cell. 8.2: Snapshots from 24 hours (gray) to 48 hours (black) taken every 30 
minutes of the profile of the insoluble pool over the cell. 8.3: The distribution of the 
keratin material between the soluble and insoluble pools over time. 8.4: The profiles of 
assembly and disassembly maximal rates used with the Michaelis-Menten type turnover 
terms in Scenario 21. The estimates are kass = 9.3819/xM/5, ks = 570.73//M, 
kdis = 0.9998//M/s leading to kmax = 1.976/iM/s in ( |l6| ) and kj = 976.07//M/s. 
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Figure 9. Estimate of the space-dependent speed. Estimate of the 


space-dependent speed given in (13) for the variable speed u(x) of the active transport 
of the insoluble pool. 
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Figure 10. Profile for localized assembly rate of type Sources. The profile for 
kassix) defined in (141 with kmax = I, obtained by fitting the profile of assembly 

and shown in Fig. [^3. 


regions Sources published in 


21 



-20 -10 0 10 20 


Location (p m) 

Figure 11. Profile for localized disassembly rate of type Sinks. The profile for 

12 ; = 1, obtained by fitting the profile of disassembly 
and shown in Fig. [^3. 
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